/** \file    diff.h
  * \brief   Simplest differentiation approximations for rasters.
  * \date    2002
  * \author  Fedor Chelnokov
  * \version $Header: /home/cscom1r4/cvsroot/cs/include/raster/diff.h,v 1.1 2004/03/26 05:58:02 cscom1r4 Exp $
  */

#ifndef _diff_h_
#define _diff_h_

namespace raster {

//////////////////////////////////////////////////////////////
//smth			   represents appropriate derivative for point
//dx(x,y)        - (x+0.5,y)
//dx(dx(x,y))    - (x+1.0,y)
//dx(dx(dx(x,y)))- (x+1.5,y)
//////////////////////////////////////////////////////////////

template <class I, class O>
void dx(const I& i, O &o)
{
	assert(i.width() == o.width() && i.height() == o.height());
	size_t w1 = i.width() - 1;
	I::const_iterator ii;
	for (ii = i.begin(); ii != i.end(); ++ii) {
		if (ii.x() == w1)
			o[ii] = O::elem_type();
		else
			o[ii] = i(ii.x()+1, ii.y()) - *ii; 
	}
}

template <class I, class O>
void dy(const I& i, O &o)
{
	assert(i.width() == o.width() && i.height() == o.height());
	size_t h1 = i.height() - 1;
	I::const_iterator ii;
	for (ii = i.begin(); ii != i.end(); ++ii) {
		if (ii.y() == h1)
			o[ii] = O::elem_type();
		else
			o[ii] = i(ii.x(), ii.y()+1) - *ii; 
	}
}

//////////////////////////////////////////////////////////////
//md?(x,y) represents first derivative for point (x+0.5,y+0.5)
//////////////////////////////////////////////////////////////

template <class I, class O>
void mdx(const I& i, O &o)
{
	assert(i.width() == o.width() && i.height() == o.height());
	size_t w1 = i.width() - 1, h1 = i.height() - 1;
	I::const_iterator ii;
	for (ii = i.begin(); ii != i.end(); ++ii) {
		size_t x = ii.x(), y = ii.y();
		if (x == w1 || y == h1)
			o[ii] = O::elem_type();
		else
			o[ii] = i(x+1,y) - i(x,y) + i(x+1,y+1) - i(x,y+1); 
	}
}

template <class I, class O>
void mdy(const I& i, O &o)
{
	assert(i.width() == o.width() && i.height() == o.height());
	size_t w1 = i.width() - 1, h1 = i.height() - 1;
	I::const_iterator ii;
	for (ii = i.begin(); ii != i.end(); ++ii) {
		size_t x = ii.x(), y = ii.y();
		if (x == w1 || y == h1)
			o[ii] = O::elem_type();
		else
			o[ii] = i(x,y+1) - i(x,y) + i(x+1,y+1) - i(x+1,y); 
	}
}

//////////////////////////////////////////////////////////////
//symmetric approximations of second order
//sd?(x,y) represents first derivative for point (x,y)
//////////////////////////////////////////////////////////////

template <class I, class O>
void sdx(const I& i, O &o)
{
	assert(i.width() == o.width() && i.height() == o.height());
	size_t w1 = i.width() - 1;
	I::const_iterator ii;
	for (ii = i.begin(); ii != i.end(); ++ii) {
		if (ii.x() == 0 || ii.x() == w1)
			o[ii] = O::elem_type();
		else
			o[ii] = i(ii.x()+1, ii.y()) - i(ii.x()-1, ii.y()); 
	}
}

template <class I, class O>
void sdy(const I& i, O &o)
{
	assert(i.width() == o.width() && i.height() == o.height());
	size_t h1 = i.height() - 1;
	I::const_iterator ii;
	for (ii = i.begin(); ii != i.end(); ++ii) {
		if (ii.y() == 0 || ii.y() == h1)
			o[ii] = O::elem_type();
		else
			o[ii] = i(ii.x(), ii.y()+1) - i(ii.x(), ii.y()-1); 
	}
}

} //namespace raster

#endif //_diff_h_
